function [x, u, t, info] = burgers(varargin)

%  [x,u,t,info] = burgers;


%BURGERS Solve the viscous Burgers equation by an implicit finite difference method.
%
%   [X,U,T,INFO] = BURGERS() solves
%
%       u_t + u*u_x = epsilon*u_xx,       xL < x < xR,
%
%   subject to
%
%       u(xL,t) = a0,    u(xR,t) = b0,
%
%   and the step initial condition
%
%       u(x,0) = a0,  x < xStep,
%              = b0,  x >= xStep.
%
%   Backward Euler is used in time, centered differences are used in
%   space, and the nonlinear equations at each time step are solved by
%   Newton's method. Each Newton correction is computed with a Thomas
%   tridiagonal solve. This is a MATLAB translation of burgers.f.
%
%   The default parameters reproduce the Fortran program:
%
%       xL = 0, xR = 1, a0 = 1, b0 = -1, xStep = 0.25,
%       nx = 10000, nt = 10000, tmax = 10000, epsilon = 0.01.
%
%   Name-value options are:
%
%       'xL', 'xR', 'a0', 'b0', 'xStep', 'nx', 'nt', 'tmax',
%       'epsilon', 'Tolerance', 'MaxNewtonIterations',
%       'OutputFile', and 'MakePlot'.
%
%   Example using a smaller test calculation:
%
%       [x,u,t] = burgers('nx',400,'nt',200,'tmax',1,'MakePlot',true);
%
%   Calling BURGERS without output arguments plots the final solution.
%
%   Note: The original values nx = nt = 10000 require a substantial
%   amount of computation in MATLAB.

    p = inputParser;
    p.FunctionName = 'burgers';

    addParameter(p,'xL',0,@isFiniteScalar);
    addParameter(p,'xR',1,@isFiniteScalar);
    addParameter(p,'a0',1,@isFiniteScalar);
    addParameter(p,'b0',-1,@isFiniteScalar);
    addParameter(p,'xStep',0.25,@isFiniteScalar);
    addParameter(p,'nx',10000,@isPositiveIntegerScalar);
    addParameter(p,'nt',10000,@isPositiveIntegerScalar);
    addParameter(p,'tmax',10000,@isNonnegativeFiniteScalar);
    addParameter(p,'epsilon',0.01,@isPositiveFiniteScalar);
    addParameter(p,'Tolerance',1.0e-6,@isPositiveFiniteScalar);
    addParameter(p,'MaxNewtonIterations',100,@isPositiveIntegerScalar);
    addParameter(p,'OutputFile','data_10000.txt',@isTextScalar);
    addParameter(p,'MakePlot',false,@isLogicalScalar);

    parse(p,varargin{:});
    prm = p.Results;

    if prm.xR <= prm.xL
        error('burgers:InvalidInterval','xR must be greater than xL.');
    end
    if prm.xStep < prm.xL || prm.xStep > prm.xR
        error('burgers:InvalidStepLocation', ...
              'xStep must lie in the interval [xL,xR].');
    end
    if prm.nt < 2
        error('burgers:InvalidTimeGrid','nt must be at least 2.');
    end

    nx = prm.nx;
    nt = prm.nt;
    ep = prm.epsilon;

    dx = (prm.xR - prm.xL)/(nx + 1);
    dx2 = dx^2;
    dt = prm.tmax/(nt - 1);

    if dt == 0
        error('burgers:ZeroTimeStep', ...
              'tmax must be positive when nt is greater than one.');
    end

    % The vector includes both boundary points:
    % x(1) = xL and x(nx+2) = xR.
    x = prm.xL + (0:nx+1)'*dx;

    % Step initial condition from the Fortran function phi.
    u = prm.b0*ones(nx + 2,1);
    u(x < prm.xStep) = prm.a0;
    u(1) = prm.a0;
    u(end) = prm.b0;

    totalNewtonIterations = 0;
    maximumNewtonIterations = 0;

    for it = 2:nt
        oldU = u;

        % In the Fortran program, q = -oldU/(dt*epsilon).
        q = -oldU(2:nx+1)/(dt*ep);

        % Use the preceding time level as the initial Newton iterate.
        y = oldU;
        y(1) = prm.a0;
        y(end) = prm.b0;

        converged = false;

        for newtonIteration = 1:prm.MaxNewtonIterations
            center = y(2:nx+1);
            left   = y(1:nx);
            right  = y(3:nx+2);

            z = (right - left)/(2*dx);

            % Equation written as y_xx = f(x,y,y_x).
            f  = (center.*z + center/dt)/ep + q;
            fy = (z + 1/dt)/ep;
            fz = center/ep;

            residual = 2*center - right - left + dx2*f;

            diagonal = 2 + dx2*fy;
            lowerRow = -1 - 0.5*dx*fz;
            upperRow = -1 + 0.5*dx*fz;

            correction = tridiagonalSolve( ...
                lowerRow(2:end), diagonal, upperRow(1:end-1), -residual);

            y(2:nx+1) = center + correction;

            if norm(correction,inf) <= prm.Tolerance
                converged = true;
                break
            end
        end

        totalNewtonIterations = totalNewtonIterations + newtonIteration;
        maximumNewtonIterations = max(maximumNewtonIterations, ...
                                      newtonIteration);

        if ~converged
            error('burgers:NewtonNotConverged', ...
                 ['Newton iteration failed at time step %d (t = %.16g) ' ...
                  'after %d iterations.'], ...
                  it, (it - 1)*dt, prm.MaxNewtonIterations);
        end

        u = y;
    end

    t = prm.tmax;

    outputFile = char(prm.OutputFile);
    if ~isempty(outputFile)
        fid = fopen(outputFile,'w');
        if fid == -1
            error('burgers:FileOpenError', ...
                  'Unable to open "%s" for writing.',outputFile);
        end
        cleanupObject = onCleanup(@() fclose(fid)); 
        fprintf(fid,'%14.5e%14.5e\n',[x,u].');
    end

    info = struct( ...
        'dx',dx, ...
        'dt',dt, ...
        'epsilon',ep, ...
        'timeSteps',nt - 1, ...
        'totalNewtonIterations',totalNewtonIterations, ...
        'maximumNewtonIterations',maximumNewtonIterations, ...
        'outputFile',outputFile);

    if prm.MakePlot || nargout == 0
        plot(x,u,'LineWidth',1.25)
        xlabel('x')
        ylabel('u(x,t)')
        title(sprintf('Viscous Burgers equation at t = %g',t))
        grid on
    end
end

function solution = tridiagonalSolve(lower,diagonal,upper,rhs)
%TRIDIAGONALSOLVE Solve a tridiagonal linear system by Thomas elimination.

    n = numel(diagonal);

    if numel(rhs) ~= n || numel(lower) ~= max(n - 1,0) || ...
            numel(upper) ~= max(n - 1,0)
        error('burgers:InvalidTridiagonalDimensions', ...
              'The tridiagonal arrays have inconsistent dimensions.');
    end

    modifiedUpper = zeros(max(n - 1,0),1);
    modifiedRhs = zeros(n,1);

    pivotScale = max(1,norm(diagonal,inf));
    pivotTolerance = eps(pivotScale);

    pivot = diagonal(1);
    checkPivot(pivot,pivotTolerance);

    if n > 1
        modifiedUpper(1) = upper(1)/pivot;
    end
    modifiedRhs(1) = rhs(1)/pivot;

    for j = 2:n
        pivot = diagonal(j) - lower(j - 1)*modifiedUpper(j - 1);
        checkPivot(pivot,pivotTolerance);

        if j < n
            modifiedUpper(j) = upper(j)/pivot;
        end
        modifiedRhs(j) = ...
            (rhs(j) - lower(j - 1)*modifiedRhs(j - 1))/pivot;
    end

    solution = zeros(n,1);
    solution(n) = modifiedRhs(n);

    for j = n-1:-1:1
        solution(j) = modifiedRhs(j) - modifiedUpper(j)*solution(j + 1);
    end
end

function checkPivot(pivot,pivotTolerance)
%CHECKPIVOT Guard against division by a zero or numerically tiny pivot.

    if ~isfinite(pivot) || abs(pivot) <= pivotTolerance
        error('burgers:ZeroPivot', ...
              'A zero or numerically singular pivot occurred.');
    end
end

function tf = isFiniteScalar(value)
    tf = isnumeric(value) && isreal(value) && isscalar(value) && ...
         isfinite(value);
end

function tf = isNonnegativeFiniteScalar(value)
    tf = isFiniteScalar(value) && value >= 0;
end

function tf = isPositiveFiniteScalar(value)
    tf = isFiniteScalar(value) && value > 0;
end

function tf = isPositiveIntegerScalar(value)
    tf = isFiniteScalar(value) && value >= 1 && value == floor(value);
end

function tf = isTextScalar(value)
    tf = ischar(value) || (isstring(value) && isscalar(value));
end

function tf = isLogicalScalar(value)
    tf = (islogical(value) || isnumeric(value)) && isscalar(value) && ...
         isfinite(value) && (value == 0 || value == 1);
end
